source("~/Dropbox/R_package_testing/brmsish/R/extended_summary.R")
source("~/Dropbox/R_package_testing/brmsish/R/helpers.R")
label_spectro_temp <- function(wave, reference = NULL, detection = NULL, envelope = FALSE,
threshold = NULL, smooth = 5, collevels = seq(-100, 0, 5), palette = viridis::viridis,
template.correlation = NULL, line.x.position = 2, hop.size = NULL, col.line = NULL,
...) {
# adjust wl based on hope.size
if (!is.null(hop.size))
wl <- round(wave@samp.rate * hop.size/1000, 0)
# reset graphic device on exit
oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
if (envelope | !is.null(template.correlation))
par(mfrow = c(2, 1), mar = c(0, 4, 1, 1)) else par(mar = c(4, 4, 1, 1))
# plot spectrogram
seewave::spectro(wave = wave, grid = FALSE, scale = FALSE, palette = palette,
collevels = collevels, axisX = if (envelope | !is.null(template.correlation))
FALSE else TRUE, ...)
# plot detection
if (!is.null(reference))
for (i in seq_len(nrow(reference))) lines(x = (reference[i, c("start", "end")]),
y = rep(line.x.position, 2), col = col.line[1], lwd = 7, lend = 2)
# plot detection
if (!is.null(detection))
for (i in seq_len(nrow(detection))) lines(x = (detection[i, c("start", "end")]),
y = rep(line.x.position - 0.3, 2), col = col.line[2], lwd = 7, lend = 2)
usr <- par("usr")
# add legend
if (!is.null(detection) & !is.null(reference))
legend(x = usr[2] * 0.98, y = usr[4] * 0.98, col = col.line[1:2], legend = c("Referencia",
"Detección"), text.width = 0.3, lwd = 4, bg = "#FFFFFFE6", xjust = 1,
yjust = 1)
if (is.null(detection) & !is.null(reference))
legend(x = usr[2] * 0.98, y = usr[4] * 0.98, col = col.line[1], text.width = 0.3,
legend = c("Referencia"), lwd = 4, bg = "#FFFFFFE6", xjust = 1, yjust = 1)
if (is.null(reference) & !is.null(detection))
legend(x = usr[2] * 0.98, y = usr[4] * 0.98, col = col.line[2], legend = c("Detección"),
lwd = 4, text.width = 0.3, bg = "#FFFFFFE6", xjust = 1, yjust = 1)
if (envelope) {
# set graphic device for envelope
par(mar = c(4, 4, 0.3, 1))
if (!is.null(smooth))
smooth <- round(wave@samp.rate * smooth/1000, 0)
# plot envelope
seewave::env(wave, colwave = "#07889B", ssmooth = smooth)
# add threshold line
if (!is.null(threshold))
abline(h = par("usr")[4] * threshold/100, col = col.line[3], lwd = 3)
} else if (!is.null(template.correlation)) {
# set graphic device for correlations
par(mar = c(4, 4, 0.3, 1))
plot(x = seq(template.correlation$template.duration/2, duration(wave) - template.correlation$template.duration/2,
length.out = length(template.correlation$correlation.scores)), y = template.correlation$correlation.scores,
type = "l", xlab = "Tiempo (s)", ylab = "Correlación", col = "#07889B",
lwd = 1.6, xaxs = "i", xlim = c(0, duration(wave)))
# add threshold line
if (!is.null(threshold))
abline(h = threshold, col = col.line[3], lwd = 3)
}
}
clim_dat_2020 <- read_excel("./data/datos_metereologicos/Estacion_met/Dat_met_estacion_2022-01-11.xlsx",
sheet = "2020")
clim_dat_2019 <- read_excel("./data/datos_metereologicos/Estacion_met/Dat_met_estacion_2022-01-11.xlsx",
sheet = "2019")
clim_dat <- rbind(clim_dat_2019, clim_dat_2020)
clim_dat <- clim_dat[, c("filename", "Año", "Mes", "Día", "Hora", "Temp (°C)",
"Humedad Relat.", "Precipitación")]
names(clim_dat) <- c("filename", "year", "month", "day", "hour", "temp", "HR", "rain")
clim_dat <- aggregate(cbind(rain, temp, HR) ~ filename + year + month + day + hour,
clim_dat, mean)
clim_dat$year <- clim_dat$year + 2000
clim_dat$date <- as.Date(paste(clim_dat$year, clim_dat$month, clim_dat$day, sep = "-"))
clim_dat$date_hour <- paste(gsub("-", "", clim_dat$date), clim_dat$hour, sep = "-")
call_dat <- read.csv("./data/processed/call_rate_per_date_time_and_site.csv")
# remove 5 pm data and keep only form Sukia
call_dat <- call_dat[call_dat$hour != 17, ]
call_dat <- call_dat[call_dat$site != "LAGCHIMU", ]
call_dat$date_hour <- paste(sapply(as.character(call_dat$site_date_hour), function(x) strsplit(x,
"_")[[1]][2]), call_dat$hour, sep = "-")
call_dat$temp <- sapply(1:nrow(call_dat), function(x) {
y <- clim_dat$temp[clim_dat$date_hour == call_dat$date_hour[x]]
if (length(y) < 1)
y <- NA
return(y)
})
call_dat$HR <- sapply(1:nrow(call_dat), function(x) {
y <- clim_dat$HR[clim_dat$date_hour == call_dat$date_hour[x]]
if (length(y) < 1)
y <- NA
return(y)
})
call_dat$rain <- sapply(1:nrow(call_dat), function(x) {
y <- clim_dat$rain[clim_dat$date_hour == call_dat$date_hour[x]]
if (length(y) < 1)
y <- NA
return(y)
})
# proportion of acoustic data with climatic data
sum(call_dat$date_hour %in% clim_dat$date_hour)/nrow(call_dat)
sum(!is.na(call_dat$temp))/nrow(call_dat)
call_dat <- call_dat[!is.na(call_dat$temp), ]
call_dat$day <- as.numeric(substr(sapply(as.character(call_dat$site_date_hour), function(x) strsplit(x,
"_")[[1]][2]), 7, 8))
call_dat$date <- as.Date(paste(call_dat$year, call_dat$month, call_dat$day, sep = "-"))
call_dat$moon.date <- ifelse(call_dat$hour < 12, as.Date(call_dat$date - 1), as.Date(call_dat$date))
call_dat$moon.date <- as.Date(call_dat$moon.date, origin = "1970-01-02")
## add moon
call_dat$moonlight <- lunar.illumination(call_dat$moon.date, shift = -6)
call_dat$date_hour_min <- strptime(paste(paste(call_dat$year, call_dat$month, call_dat$day,
sep = "-"), paste(call_dat$hour, "00", sep = ":")), format = "%Y-%m-%d %H:%M")
call_dat$hour_diff <- as.numeric(call_dat$date_hour_min - min(call_dat$date_hour_min))/3600
call_dat$rain_24 <- sapply(1:nrow(call_dat), function(x) sum(clim_dat$rain[strptime(clim_dat$date,
format = "%Y-%m-%d") == (strptime(call_dat$date[x], format = "%Y-%m-%d") - 60 *
60 * 24)]))
call_dat$rain_48 <- sapply(1:nrow(call_dat), function(x) sum(clim_dat$rain[strptime(clim_dat$date,
format = "%Y-%m-%d") == (strptime(call_dat$date[x], format = "%Y-%m-%d") - 60 *
60 * 48)]))
clim_dat$date_hour_min <- strptime(paste(paste(clim_dat$year, clim_dat$month, clim_dat$day,
sep = "-"), paste(clim_dat$hour, "00", sep = ":")), format = "%Y-%m-%d %H:%M")
clim_dat$hour_diff <- as.numeric(clim_dat$date_hour_min - min(call_dat$date_hour_min))/3600
call_dat$prev_temp <- sapply(1:nrow(call_dat), function(x) {
# if(call_dat$hour_diff[x] < 48) pt <- NA else
pt <- mean(clim_dat$temp[clim_dat$hour_diff %in% (call_dat$hour_diff[x] - 48):(call_dat$hour_diff[x] -
24)])
return(pt)
})
write.csv(call_dat, "./data/processed/acoustic_and_climatic_data_by_hour.csv")
Purpose
- Evaluate effect of enviromental factors on vocal activity of A.
lemur
Prepare data
Descriptive stats
call_dat_site <- read.csv("./data/processed/call_rate_per_date_time_and_site.csv")
- Total number of recordings: 9681
- Total recordings per site:
agg_recs <- aggregate(rec_time ~ site, data = call_dat_site, length)
names(agg_recs)[1:2] <- c("site", "rec_count")
# print table as kable
kb <- kable(agg_recs, row.names = TRUE, digits = 3)
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
print(kb)
|
|
site
|
rec_count
|
|
1
|
LAGCHIMU
|
5127
|
|
2
|
SUKIA
|
4554
|
agg_recs <- aggregate(rec_time ~ site, data = call_dat_site, sum)
names(agg_recs)[1:2] <- c("site", "recording_time")
agg_recs$recording_time <- round((agg_recs$recording_time)/60)
# print table as kable
kb <- kable(agg_recs, row.names = TRUE, digits = 3)
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
print(kb)
|
|
site
|
recording_time
|
|
1
|
LAGCHIMU
|
1707
|
|
2
|
SUKIA
|
1517
|
agg_recs <- aggregate(n_call ~ site, data = call_dat_site, sum)
names(agg_recs)[1:2] <- c("calls", "count")
# print table as kable
kb <- kable(agg_recs, row.names = TRUE, digits = 3)
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
print(kb)
|
|
calls
|
count
|
|
1
|
LAGCHIMU
|
169485
|
|
2
|
SUKIA
|
370718
|
Call rate: 18.598587
Call rate per site:
agg_recs <- aggregate(call_rate ~ site, data = call_dat_site, mean)
agg_recs$sd <- aggregate(call_rate ~ site, data = call_dat_site, sd)[, 2]
names(agg_recs)[1:3] <- c("site", "call_rate", "sd")
# print table as kable
kb <- kable(agg_recs, row.names = TRUE, digits = 3)
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
print(kb)
|
|
site
|
call_rate
|
sd
|
|
1
|
LAGCHIMU
|
11.018
|
18.039
|
|
2
|
SUKIA
|
27.133
|
87.778
|
call_rate_hour <- read.csv("./data/processed/acoustic_and_climatic_data_by_hour.csv")
agg <- aggregate(cbind(temp, prev_temp, HR, rain, rain_24, rain_48, moonlight) ~
1, call_rate_hour, function(x) round(c(mean(x), sd(x), min(x), max(x)), 3))
agg <- as.data.frame(matrix(unlist(agg), ncol = 4, byrow = TRUE, dimnames = list(c("Temperature",
"Previous temperature", "Relative humidity", "Night rain", "Rain 24 hours", "Rain 48 hours",
"Moonlight"), c("mean", "sd", "min", "max"))))
# print table as kable
kb <- kable(agg, row.names = TRUE, digits = 3)
kb <- kable_styling(kb, bootstrap_options = c("striped", "hover", "condensed", "responsive"))
print(kb)
|
|
mean
|
sd
|
min
|
max
|
|
Temperature
|
24.229
|
2.126
|
19.472
|
31.763
|
|
Previous temperature
|
24.097
|
0.988
|
20.005
|
27.111
|
|
Relative humidity
|
90.056
|
8.464
|
55.158
|
99.940
|
|
Night rain
|
0.079
|
0.457
|
0.000
|
10.417
|
|
Rain 24 hours
|
1.441
|
3.161
|
0.000
|
25.620
|
|
Rain 48 hours
|
1.453
|
3.170
|
0.000
|
25.620
|
|
Moonlight
|
0.495
|
0.355
|
0.000
|
1.000
|
Mean and sd temperature: 24.229179
Mean previous temperature: 24.097073
Mean temperature: 24.229179
Mean cumulative rain per hour: 0.079232
Mean cumulative rain per hour previous 24 hours: 1.44064
Mean daily cumulative rain per hour previous 48 hours: 1.45317
Bayesian regression models
Scale variables and set model parameters
call_rate_hour <- read.csv("./data/processed/acoustic_and_climatic_data_by_hour.csv")
# make hour a factor
call_rate_hour$hour <- factor(call_rate_hour$hour)
# scale and mean-center
call_rate_hour$sc_temp <- scale(call_rate_hour$temp)
call_rate_hour$sc_HR <- scale(call_rate_hour$HR)
call_rate_hour$sc_rain <- scale(call_rate_hour$rain)
call_rate_hour$sc_rain_24 <- scale(call_rate_hour$rain_24)
call_rate_hour$sc_rain_48 <- scale(call_rate_hour$rain_48)
call_rate_hour$sc_moonlight <- scale(call_rate_hour$moonlight)
call_rate_hour$sc_prev_temp <- scale(call_rate_hour$prev_temp)
priors <- c(prior(normal(0, 4), class = "b"))
chains <- 4
iter <- 10000
Global models
The increase in relative humidity, decrease in temperature, increase
in the previous accumulated rain, decrease in the night rain, decrease
in the percentage of the moon illuminated cause the activity of A. lemur
to increase?
fit_glob1 <- brm(n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_moonlight +
sc_rain + sc_rain_24 + sc_prev_temp + ar(p = 2, time = hour_diff, gr = hour),
data = call_rate_hour, iter = iter, chains = chains, cores = chains, family = negbinomial(),
prior = priors, file = "./data/processed/regression_models/global_rain_24", file_refit = "always")
fit_glob2 <- brm(n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_moonlight +
sc_rain + sc_rain_48 + sc_prev_temp + ar(p = 2, time = hour_diff, gr = hour),
data = call_rate_hour, iter = iter, chains = chains, cores = chains, family = negbinomial(),
prior = priors, file = "./data/processed/regression_models/global_rain_48", file_refit = "always")
extended_summary(read.file = "./data/processed/regression_models/global_rain_24.rds",
n.posterior = 1000)
global_rain_24
|
|
priors
|
formula
|
iterations
|
chains
|
thinning
|
warmup
|
diverg_transitions
|
rhats > 1.05
|
min_bulk_ESS
|
min_tail_ESS
|
seed
|
|
1
|
b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0,
2.5) shape-gamma(0.01, 0.01)
|
n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_moonlight + sc_rain
+ sc_rain_24 + sc_prev_temp + ar(p = 2, time = hour_diff, gr = hour)
|
10000
|
4
|
1
|
5000
|
0
|
0
|
9377.5
|
12503.8
|
536473866
|
|
|
Estimate
|
l-95% CI
|
u-95% CI
|
Rhat
|
Bulk_ESS
|
Tail_ESS
|
|
b_Intercept
|
0.698
|
0.636
|
0.760
|
1
|
9377.5
|
12503.8
|
|
b_sc_temp
|
0.573
|
0.493
|
0.654
|
1
|
10930.3
|
12951.7
|
|
b_sc_HR
|
0.307
|
0.233
|
0.381
|
1
|
14333.7
|
14034.8
|
|
b_sc_moonlight
|
-0.090
|
-0.149
|
-0.031
|
1
|
12403.9
|
13102.6
|
|
b_sc_rain
|
0.038
|
0.002
|
0.075
|
1
|
18321.1
|
15080.0
|
|
b_sc_rain_24
|
0.091
|
0.053
|
0.130
|
1
|
18798.2
|
14951.0
|
|
b_sc_prev_temp
|
-0.024
|
-0.069
|
0.022
|
1
|
14546.5
|
13679.8
|

extended_summary(read.file = "./data/processed/regression_models/global_rain_48.rds",
n.posterior = 1000)
global_rain_48
|
|
priors
|
formula
|
iterations
|
chains
|
thinning
|
warmup
|
diverg_transitions
|
rhats > 1.05
|
min_bulk_ESS
|
min_tail_ESS
|
seed
|
|
1
|
b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0,
2.5) shape-gamma(0.01, 0.01)
|
n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_moonlight + sc_rain
+ sc_rain_48 + sc_prev_temp + ar(p = 2, time = hour_diff, gr = hour)
|
10000
|
4
|
1
|
5000
|
0
|
0
|
12375.1
|
13473.4
|
58053599
|
|
|
Estimate
|
l-95% CI
|
u-95% CI
|
Rhat
|
Bulk_ESS
|
Tail_ESS
|
|
b_Intercept
|
0.703
|
0.640
|
0.766
|
1
|
12375.1
|
13699.9
|
|
b_sc_temp
|
0.566
|
0.485
|
0.649
|
1
|
13593.6
|
13473.4
|
|
b_sc_HR
|
0.303
|
0.228
|
0.379
|
1
|
17757.9
|
15101.0
|
|
b_sc_moonlight
|
-0.098
|
-0.158
|
-0.039
|
1
|
18857.2
|
14680.5
|
|
b_sc_rain
|
0.031
|
-0.006
|
0.068
|
1
|
22865.9
|
15592.2
|
|
b_sc_rain_48
|
-0.006
|
-0.044
|
0.033
|
1
|
24427.4
|
16403.0
|
|
b_sc_prev_temp
|
-0.037
|
-0.084
|
0.010
|
1
|
20592.6
|
15806.9
|

Temperature
The increase in temperature at night causes that activity of A. lemur
to decrease?
fit2.1 <- brm(n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_24 + ar(p = 2,
time = hour_diff, gr = hour), data = call_rate_hour, iter = iter, chains = chains,
cores = chains, family = negbinomial(), prior = priors, file = "./data/processed/regression_models/temp2.1",
file_refit = "always")
fit2.2 <- brm(n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_48 + ar(p = 2,
time = hour_diff, gr = hour), data = call_rate_hour, iter = iter, chains = chains,
cores = chains, family = negbinomial(), prior = priors, file = "./data/processed/regression_models/temp2.2",
file_refit = "always")
fit2.3 <- brm(n_call | resp_rate(rec_time) ~ sc_prev_temp + ar(p = 2, time = hour_diff,
gr = hour), data = call_rate_hour, iter = iter, chains = chains, cores = chains,
family = negbinomial(), prior = priors, file = "./data/processed/regression_models/temp2.3",
file_refit = "always")
extended_summary(read.file = "./data/processed/regression_models/temp2.1.rds", n.posterior = 1000)
temp2.1
|
|
priors
|
formula
|
iterations
|
chains
|
thinning
|
warmup
|
diverg_transitions
|
rhats > 1.05
|
min_bulk_ESS
|
min_tail_ESS
|
seed
|
|
1
|
b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0,
2.5) shape-gamma(0.01, 0.01)
|
n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_24 + ar(p =
2, time = hour_diff, gr = hour)
|
10000
|
4
|
1
|
5000
|
1
|
0
|
10328.4
|
12111.8
|
1846130948
|
|
|
Estimate
|
l-95% CI
|
u-95% CI
|
Rhat
|
Bulk_ESS
|
Tail_ESS
|
|
b_Intercept
|
0.698
|
0.634
|
0.762
|
1
|
10328.4
|
12111.8
|
|
b_sc_temp
|
0.312
|
0.261
|
0.361
|
1
|
12331.5
|
14350.8
|
|
b_sc_rain
|
0.044
|
0.008
|
0.083
|
1
|
20658.2
|
15981.9
|
|
b_sc_rain_24
|
0.090
|
0.052
|
0.128
|
1
|
20542.9
|
16470.9
|

extended_summary(read.file = "./data/processed/regression_models/temp2.2.rds", n.posterior = 1000)
temp2.2
|
|
priors
|
formula
|
iterations
|
chains
|
thinning
|
warmup
|
diverg_transitions
|
rhats > 1.05
|
min_bulk_ESS
|
min_tail_ESS
|
seed
|
|
1
|
b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0,
2.5) shape-gamma(0.01, 0.01)
|
n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_48 + ar(p =
2, time = hour_diff, gr = hour)
|
10000
|
4
|
1
|
5000
|
0
|
0
|
8873.8
|
11982.2
|
2046380001
|
|
|
Estimate
|
l-95% CI
|
u-95% CI
|
Rhat
|
Bulk_ESS
|
Tail_ESS
|
|
b_Intercept
|
0.701
|
0.635
|
0.766
|
1
|
8873.8
|
11982.2
|
|
b_sc_temp
|
0.305
|
0.255
|
0.355
|
1
|
10391.0
|
13564.9
|
|
b_sc_rain
|
0.036
|
0.000
|
0.073
|
1
|
17117.4
|
15242.3
|
|
b_sc_rain_48
|
0.006
|
-0.031
|
0.044
|
1
|
17732.8
|
15884.1
|

extended_summary(read.file = "./data/processed/regression_models/temp2.3.rds", n.posterior = 1000)
temp2.3
|
|
priors
|
formula
|
iterations
|
chains
|
thinning
|
warmup
|
diverg_transitions
|
rhats > 1.05
|
min_bulk_ESS
|
min_tail_ESS
|
seed
|
|
1
|
b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0,
2.5) shape-gamma(0.01, 0.01)
|
n_call | resp_rate(rec_time) ~ sc_prev_temp + ar(p = 2, time =
hour_diff, gr = hour)
|
10000
|
4
|
1
|
5000
|
0
|
0
|
20200
|
15620.2
|
1962664000
|
|
|
Estimate
|
l-95% CI
|
u-95% CI
|
Rhat
|
Bulk_ESS
|
Tail_ESS
|
|
b_Intercept
|
0.675
|
0.605
|
0.745
|
1
|
20200.0
|
15620.2
|
|
b_sc_prev_temp
|
-0.002
|
-0.048
|
0.044
|
1
|
36045.7
|
16766.2
|

Relative humidity
Does an increase in relative humidity cause the activity of A. lemur
to increase?
fit.3.1 <- brm(n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_rain + sc_rain_24 +
ar(p = 2, time = hour_diff, gr = hour), data = call_rate_hour, iter = iter, chains = chains,
cores = chains, family = negbinomial(), prior = priors, file = "./data/processed/regression_models/RH3.1",
file_refit = "always")
fit.3.2 <- brm(n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_rain + sc_rain_48 +
ar(p = 2, time = hour_diff, gr = hour), data = call_rate_hour, iter = iter, chains = chains,
cores = chains, family = negbinomial(), prior = priors, file = "./data/processed/regression_models/RH3.2",
file_refit = "always") # Este tengo que doblemente confirmarlo
extended_summary(read.file = "./data/processed/regression_models/RH3.1.rds", n.posterior = 1000)
RH3.1
|
|
priors
|
formula
|
iterations
|
chains
|
thinning
|
warmup
|
diverg_transitions
|
rhats > 1.05
|
min_bulk_ESS
|
min_tail_ESS
|
seed
|
|
1
|
b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0,
2.5) shape-gamma(0.01, 0.01)
|
n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_rain + sc_rain_24 +
ar(p = 2, time = hour_diff, gr = hour)
|
10000
|
4
|
1
|
5000
|
0
|
0
|
10000
|
12606.5
|
1854758574
|
|
|
Estimate
|
l-95% CI
|
u-95% CI
|
Rhat
|
Bulk_ESS
|
Tail_ESS
|
|
b_Intercept
|
0.699
|
0.637
|
0.760
|
1
|
10000.0
|
12606.5
|
|
b_sc_temp
|
0.578
|
0.499
|
0.656
|
1
|
13756.3
|
13884.5
|
|
b_sc_HR
|
0.314
|
0.242
|
0.385
|
1
|
17498.1
|
14651.1
|
|
b_sc_rain
|
0.037
|
0.002
|
0.075
|
1
|
24287.1
|
16414.2
|
|
b_sc_rain_24
|
0.096
|
0.058
|
0.135
|
1
|
22809.6
|
15579.0
|

extended_summary(read.file = "./data/processed/regression_models/RH3.2.rds", n.posterior = 1000)
RH3.2
|
|
priors
|
formula
|
iterations
|
chains
|
thinning
|
warmup
|
diverg_transitions
|
rhats > 1.05
|
min_bulk_ESS
|
min_tail_ESS
|
seed
|
|
1
|
b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0,
2.5) shape-gamma(0.01, 0.01)
|
n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_rain + sc_rain_48 +
ar(p = 2, time = hour_diff, gr = hour)
|
10000
|
4
|
1
|
5000
|
0
|
0
|
8492.95
|
12550.7
|
1444793977
|
|
|
Estimate
|
l-95% CI
|
u-95% CI
|
Rhat
|
Bulk_ESS
|
Tail_ESS
|
|
b_Intercept
|
0.703
|
0.641
|
0.765
|
1
|
8492.95
|
12697.7
|
|
b_sc_temp
|
0.567
|
0.488
|
0.648
|
1
|
9626.63
|
12550.7
|
|
b_sc_HR
|
0.308
|
0.235
|
0.382
|
1
|
12079.96
|
13794.6
|
|
b_sc_rain
|
0.029
|
-0.007
|
0.066
|
1
|
18138.82
|
15107.9
|
|
b_sc_rain_48
|
0.005
|
-0.032
|
0.043
|
1
|
16926.99
|
14915.2
|

Moon
Decreasing the percentage of the moon illuminated causes an increase
in A. lemur activity?
fit.4 <- brm(n_call | resp_rate(rec_time) ~ sc_moonlight + ar(p = 2, time = hour_diff,
gr = hour), data = call_rate_hour, iter = iter, chains = chains, cores = chains,
family = negbinomial(), prior = priors, file = "./data/processed/regression_models/moon4",
file_refit = "always")
extended_summary(read.file = "./data/processed/regression_models/moon4.rds", n.posterior = 1000)
moon4
|
|
priors
|
formula
|
iterations
|
chains
|
thinning
|
warmup
|
diverg_transitions
|
rhats > 1.05
|
min_bulk_ESS
|
min_tail_ESS
|
seed
|
|
1
|
b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0,
2.5) shape-gamma(0.01, 0.01)
|
n_call | resp_rate(rec_time) ~ sc_moonlight + ar(p = 2, time =
hour_diff, gr = hour)
|
10000
|
4
|
1
|
5000
|
0
|
0
|
19357
|
14632
|
1734290538
|
|
|
Estimate
|
l-95% CI
|
u-95% CI
|
Rhat
|
Bulk_ESS
|
Tail_ESS
|
|
b_Intercept
|
0.674
|
0.606
|
0.743
|
1
|
19357.0
|
14632.0
|
|
b_sc_moonlight
|
-0.128
|
-0.194
|
-0.063
|
1
|
31443.5
|
15209.9
|

Night rain
If the night rain decreases can affect the A. lemur activity to
increase?
fit.5.1 <- brm(n_call | resp_rate(rec_time) ~ sc_prev_temp + sc_rain + sc_rain_48 +
ar(p = 2, time = hour_diff, gr = hour), data = call_rate_hour, iter = iter, chains = chains,
cores = chains, family = negbinomial(), prior = priors, file = "./data/processed/regression_models/night_rain5.1",
file_refit = "always")
fit.5.2 <- brm(n_call | resp_rate(rec_time) ~ sc_prev_temp + sc_rain + sc_rain_24 +
ar(p = 2, time = hour_diff, gr = hour), data = call_rate_hour, iter = iter, chains = chains,
cores = chains, family = negbinomial(), prior = priors, file = "./data/processed/regression_models/night_rain5.2",
file_refit = "always")
extended_summary(read.file = "./data/processed/regression_models/night_rain5.1.rds",
n.posterior = 1000)
night_rain5.1
|
|
priors
|
formula
|
iterations
|
chains
|
thinning
|
warmup
|
diverg_transitions
|
rhats > 1.05
|
min_bulk_ESS
|
min_tail_ESS
|
seed
|
|
1
|
b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0,
2.5) shape-gamma(0.01, 0.01)
|
n_call | resp_rate(rec_time) ~ sc_prev_temp + sc_rain + sc_rain_48 +
ar(p = 2, time = hour_diff, gr = hour)
|
10000
|
4
|
1
|
5000
|
0
|
0
|
15449.5
|
13446.5
|
1903747053
|
|
|
Estimate
|
l-95% CI
|
u-95% CI
|
Rhat
|
Bulk_ESS
|
Tail_ESS
|
|
b_Intercept
|
0.676
|
0.607
|
0.745
|
1.001
|
15449.5
|
13446.5
|
|
b_sc_prev_temp
|
0.003
|
-0.046
|
0.052
|
1
|
27583.8
|
14753.3
|
|
b_sc_rain
|
-0.006
|
-0.042
|
0.030
|
1
|
29045.4
|
16584.0
|
|
b_sc_rain_48
|
0.014
|
-0.027
|
0.054
|
1
|
27750.2
|
15632.5
|

extended_summary(read.file = "./data/processed/regression_models/night_rain5.2.rds",
n.posterior = 1000)
night_rain5.2
|
|
priors
|
formula
|
iterations
|
chains
|
thinning
|
warmup
|
diverg_transitions
|
rhats > 1.05
|
min_bulk_ESS
|
min_tail_ESS
|
seed
|
|
1
|
b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0,
2.5) shape-gamma(0.01, 0.01)
|
n_call | resp_rate(rec_time) ~ sc_prev_temp + sc_rain + sc_rain_24 +
ar(p = 2, time = hour_diff, gr = hour)
|
10000
|
4
|
1
|
5000
|
0
|
0
|
21570
|
14952.6
|
1800703315
|
|
|
Estimate
|
l-95% CI
|
u-95% CI
|
Rhat
|
Bulk_ESS
|
Tail_ESS
|
|
b_Intercept
|
0.673
|
0.604
|
0.742
|
1
|
21570.0
|
14952.6
|
|
b_sc_prev_temp
|
0.007
|
-0.039
|
0.054
|
1
|
39064.3
|
16543.7
|
|
b_sc_rain
|
0.001
|
-0.035
|
0.038
|
1
|
36294.6
|
16198.6
|
|
b_sc_rain_24
|
0.077
|
0.038
|
0.117
|
1
|
34404.5
|
17412.9
|

Previous Rain
Does an increase in the accumulated previous rain causes that A.
lemur activity to increase?
fit.6.1 <- brm(n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_48 + ar(p = 2,
time = hour_diff, gr = hour), data = call_rate_hour, iter = iter, chains = chains,
cores = chains, family = negbinomial(), prior = priors, file = "./data/processed/regression_models/previous_rain6.1",
file_refit = "always")
fit.6.2 <- brm(n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_24 + ar(p = 2,
time = hour_diff, gr = hour), data = call_rate_hour, iter = iter, chains = chains,
cores = chains, family = negbinomial(), prior = priors, file = "./data/processed/regression_models/previous_rain6.2",
file_refit = "always")
fit.6.3 <- brm(n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_rain_24 + ar(p = 2,
time = hour_diff, gr = hour), data = call_rate_hour, iter = iter, chains = chains,
cores = chains, family = negbinomial(), prior = priors, file = "./data/processed/regression_models/night_rain6.3",
file_refit = "always")
extended_summary(read.file = "./data/processed/regression_models/previous_rain6.1.rds",
n.posterior = 1000)
previous_rain6.1
|
|
priors
|
formula
|
iterations
|
chains
|
thinning
|
warmup
|
diverg_transitions
|
rhats > 1.05
|
min_bulk_ESS
|
min_tail_ESS
|
seed
|
|
1
|
b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0,
2.5) shape-gamma(0.01, 0.01)
|
n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_48 + ar(p =
2, time = hour_diff, gr = hour)
|
10000
|
4
|
1
|
5000
|
0
|
0
|
8446.05
|
11423.8
|
20512253
|
|
|
Estimate
|
l-95% CI
|
u-95% CI
|
Rhat
|
Bulk_ESS
|
Tail_ESS
|
|
b_Intercept
|
0.701
|
0.636
|
0.764
|
1
|
8446.05
|
11423.8
|
|
b_sc_temp
|
0.304
|
0.254
|
0.355
|
1
|
10470.72
|
13123.6
|
|
b_sc_rain
|
0.036
|
0.000
|
0.074
|
1
|
16919.56
|
15417.1
|
|
b_sc_rain_48
|
0.007
|
-0.031
|
0.045
|
1
|
17069.15
|
15935.4
|

extended_summary(read.file = "./data/processed/regression_models/previous_rain6.2.rds",
n.posterior = 1000)
previous_rain6.2
|
|
priors
|
formula
|
iterations
|
chains
|
thinning
|
warmup
|
diverg_transitions
|
rhats > 1.05
|
min_bulk_ESS
|
min_tail_ESS
|
seed
|
|
1
|
b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0,
2.5) shape-gamma(0.01, 0.01)
|
n_call | resp_rate(rec_time) ~ sc_temp + sc_rain + sc_rain_24 + ar(p =
2, time = hour_diff, gr = hour)
|
10000
|
4
|
1
|
5000
|
0
|
0
|
11468.2
|
13337.9
|
1714781937
|
|
|
Estimate
|
l-95% CI
|
u-95% CI
|
Rhat
|
Bulk_ESS
|
Tail_ESS
|
|
b_Intercept
|
0.698
|
0.633
|
0.762
|
1
|
11468.2
|
13337.9
|
|
b_sc_temp
|
0.312
|
0.262
|
0.362
|
1
|
13834.2
|
13973.2
|
|
b_sc_rain
|
0.045
|
0.008
|
0.083
|
1
|
23413.4
|
15796.1
|
|
b_sc_rain_24
|
0.090
|
0.051
|
0.129
|
1
|
22867.5
|
15487.7
|

extended_summary(read.file = "./data/processed/regression_models/night_rain6.3.rds",
n.posterior = 1000)
night_rain6.3
|
|
priors
|
formula
|
iterations
|
chains
|
thinning
|
warmup
|
diverg_transitions
|
rhats > 1.05
|
min_bulk_ESS
|
min_tail_ESS
|
seed
|
|
1
|
b-normal(0, 4) Intercept-student_t(3, 3.6, 2.5) sderr-student_t(3, 0,
2.5) shape-gamma(0.01, 0.01)
|
n_call | resp_rate(rec_time) ~ sc_temp + sc_HR + sc_rain_24 + ar(p = 2,
time = hour_diff, gr = hour)
|
10000
|
4
|
1
|
5000
|
2
|
0
|
7697.31
|
11401.3
|
2086406674
|
|
|
Estimate
|
l-95% CI
|
u-95% CI
|
Rhat
|
Bulk_ESS
|
Tail_ESS
|
|
b_Intercept
|
0.699
|
0.636
|
0.760
|
1
|
7697.31
|
11401.3
|
|
b_sc_temp
|
0.572
|
0.493
|
0.650
|
1
|
8657.91
|
12183.9
|
|
b_sc_HR
|
0.318
|
0.244
|
0.391
|
1
|
10071.18
|
12501.7
|
|
b_sc_rain_24
|
0.093
|
0.055
|
0.131
|
1
|
14992.10
|
14450.6
|

Causal model combined effect size plot
Takes the effect sizes (and posteriors) from the right causal
models
coef_table <- read.csv("./data/processed/summary_causal_model_table_13-01-23.csv",
sep = "\t")
coef_table$variable <- coef_table$Label
coef_table$value <- coef_table$Estimate
coef_table$significance <- ifelse(coef_table$CI_low * coef_table$CI_high > 0, "sig",
"no.sig")
posteriors_l <- lapply(1:nrow(coef_table), function(x) {
# print(x)
X <- readRDS(coef_table$Model[x])
xdrws <- brms::as_draws(X)
post <- xdrws$`1`[[paste0("b_", coef_table$Variable[x])]]
out <- data.frame(variable = coef_table$Label[x], value = post, significance = coef_table$significance[x])
return(out)
})
posteriors <- do.call(rbind, posteriors_l)
coef_table$variable <- factor(coef_table$Label, levels = sort(unique(coef_table$Label),
FALSE))
posteriors$variable <- factor(posteriors$variable, levels = sort(unique(posteriors$variable),
TRUE))
fill_values <- c("#FFB9DF", "#FF7598")
fill_values <- adjustcolor(fill_values, alpha.f = 0.5)
# creat plots gg_dists <-
ggplot2::ggplot(data = posteriors, ggplot2::aes(y = variable, x = value, fill = significance)) +
ggplot2::geom_vline(xintercept = 0, col = "black", lty = 2) + ggdist::stat_halfeye(ggplot2::aes(x = value),
.width = c(0.95), normalize = "panels", color = "transparent") + ggplot2::scale_fill_manual(values = fill_values,
guide = "none") + ggplot2::geom_point(data = coef_table) + ggplot2::geom_errorbar(data = coef_table,
ggplot2::aes(xmin = CI_low, xmax = CI_high), width = 0) + ggplot2::scale_color_manual(values = coef) +
ggplot2::facet_wrap(~variable, scales = "free_y", ncol = 1) + ggplot2::theme_classic() +
ggplot2::theme(axis.ticks.length = ggplot2::unit(0, "pt"), plot.margin = ggplot2::margin(0,
0, 0, 0, "pt"), legend.position = "none", strip.background = ggplot2::element_blank(),
strip.text = ggplot2::element_blank()) + ggplot2::labs(x = "Effect size",
y = "Parameter") #+

# ggplot2::xlim(range(c(posteriors_by_chain$value, 0)) * plot.area.prop)
ggsave(filename = "./figures/summary_effect_sizes_pooled_from_multiple_causal_models.jpeg")
Conditional plots
Measured based on the global model
Rain and temperature
glob_mod_24 <- readRDS("./data/processed/regression_models/global_rain_24.rds")
conditions <- data.frame(sc_temp = c(`Low temperature` = -1, `Mean temperature` = 0,
`High temperature` = 1))
rain_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_rain", conditions = conditions),
plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + ylim(c(0, 520)) +
ggtitle("Current rain") + labs(x = "Rain", y = "Call activity (calls/hour)")
rain24_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_rain_24", conditions = conditions),
plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + ylim(c(0, 520)) +
ggtitle("Previous 24h rain") + labs(x = "Rain", y = "Call activity (calls/hour)")
glob_mod_48 <- readRDS("./data/processed/regression_models/global_rain_48.rds")
rain48_gg <- plot(conditional_effects(glob_mod_48, effects = "sc_rain_48", conditions = conditions),
plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + ylim(c(0, 520)) +
ggtitle("Previous 48h rain") + labs(x = "Rain", y = "Call activity (calls/hour)")
cowplot::plot_grid(rain_gg, rain24_gg, rain48_gg, nrow = 3)

Relative humidity and temperature
conditions <- data.frame(sc_temp = c(`Low temperature` = -1, `Mean temperature` = 0,
`High temperature` = 1))
hr_temp_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_HR", conditions = conditions),
plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + labs(x = "Relative humidity",
y = "Call activity (calls/hour)")
hr_temp_gg

conditions <- data.frame(sc_HR = c(`Low humidity` = -1, `Mean humidity` = 0, `High humidity` = 1))
temp_hr_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_temp", conditions = conditions),
plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + labs(x = "Temperature",
y = "Call activity (calls/hour)")
temp_hr_gg

Relative humidity and rain
conditions <- data.frame(sc_rain = c(`Low rain` = -1, `Mean rain` = 0, `High rain` = 1))
rain_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_HR", conditions = conditions),
plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + ylim(c(0, 75)) +
ggtitle("Current rain")
rain24_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_HR", conditions = conditions),
plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + ylim(c(0, 75)) +
ggtitle("Previous 24h rain")
rain48_gg <- plot(conditional_effects(glob_mod_48, effects = "sc_HR", conditions = conditions),
plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + ylim(c(0, 75)) +
ggtitle("Previous 48h rain")
cowplot::plot_grid(rain_gg, rain24_gg, rain48_gg, nrow = 3)

conditions <- data.frame(sc_HR = c(`Low humidity` = -1, `Mean humidity` = 0, `High humidity` = 1))
rain_hr_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_rain", conditions = conditions),
plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + ggtitle("Current rain") +
labs(x = "Rain", y = "Call activity (calls/hour)")
rain24_hr_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_rain_24", conditions = conditions),
plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + ggtitle("Previous 24h rain") +
labs(x = "Rain previous", y = "Call activity (calls/hour)")
rain48_hr_gg <- plot(conditional_effects(glob_mod_48, effects = "sc_rain_48", conditions = conditions),
plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + ggtitle("Previous 48h rain") +
labs(x = "Rain previous", y = "Call activity (calls/hour)")
cowplot::plot_grid(rain_hr_gg, rain24_hr_gg, rain48_hr_gg, nrow = 3)

Moonlight and temperature
conditions <- data.frame(sc_temp = c(`Low temperature` = -1, `Mean temperature` = 0,
`High temperature` = 1))
moon_temp_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_moonlight", conditions = conditions),
plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + labs(x = "Moonlight",
y = "Call activity (calls/hour)")
moon_temp_gg

conditions <- data.frame(sc_moonlight = c(`Low moonlight` = -1, `Mean moonlight` = 0,
`High moonlight` = 1))
temp_moon_gg <- plot(conditional_effects(glob_mod_24, effects = "sc_temp", conditions = conditions),
plot = FALSE)[[1]] + scale_color_viridis_d() + theme_classic() + labs(x = "Temperature",
y = "Call activity (calls/hour)")
temp_moon_gg

# st <- data.frame(sound.files = 'FIGURA_LEMUR.wav', selec = 1, start = 0.1,
# end = 0.4)
# tweak_spectro(st, length.out = 20, ovlp = 99, wl = c(100, 1000), pal =
# c('reverse.gray.colors.2', 'viridis', 'reverse.terrain.colors'), path =
# './figures/Figura_canto_lemur', flim = c(1, 4), ncol = 10, nrow = 6, width =
# 15, collev.min = c(-110))
wav <- readWave("./figures/Figura_canto_lemur/FIGURA_LEMUR.wav")
graphics.off()
par(mfrow = c(2, 1), mar = c(0, 4, 4, 1))
warbleR:::spectro_wrblr_int2(wav, flim = c(1, 4), wl = 700, grid = F, collevels = seq(-125,
0, 1), ovlp = 0, palette = viridis, axisX = FALSE)
peaks <- c(0.320165, 0.9, 1.7, 2.6, 3.5, 4.2, 4.5) # + seq(-0.1, 0.2, length.out = 7)
valleys <- seq(0, duration(wav) + 0.3, length.out = 100)
cor_dat <- data.frame(time = c(peaks, valleys), cor = c(rep(0.7, length(peaks)),
rep(0.1, length(valleys))))
cor_dat$cor <- cor_dat$cor + rnorm(nrow(cor_dat), sd = 0.03)
cor_dat$cor <- smoothw(cor_dat$cor, wl = 2, f = 10)
cor_dat <- cor_dat[order(cor_dat$time), ]
par(mar = c(2, 4, 0, 1))
plot(cor_dat$time, cor_dat$cor, type = "l", xaxs = "i")
spectro(wav, flim = c(1, 4), wl = 700, grid = F, collevels = seq(-125, 0, 1), ovlp = 99,
palette = viridis, axisX = FALSE, tlim = c(1.53, 1.63), scale = FALSE, flab = "Frecuencia (kHz)",
tlab = "Tiempo (s)")
template <- data.frame(sound.files = "FIGURA_LEMUR.wav", selec = 1, start = 1.55,
end = 1.61, bottom.freq = 2, top.freq = 3)
# get correlations
correlations <- template_correlator(templates = template, files = "FIGURA_LEMUR.wav",
path = "./figures/Figura_canto_lemur/")
thresh <- 0.7
# run detection
detection <- template_detector(template.correlations = correlations, threshold = thresh)
reference <- template_detector(template.correlations = correlations, threshold = 0.55)
detection
# plot spectrogram
label_spectro_temp(wave = wav, reference = reference, detection = detection, template.correlation = correlations[[1]],
flim = c(1, 4), threshold = thresh, hop.size = 10, ovlp = 50, collevels = seq(-125,
0, 1), col.line = c("#AFBF35", "#F20505", "#F2622E"), flab = "Frecuencia (kHz)",
tlab = "Tiempo (s)")
# 012623 034941 AFBF35 F2622E F20505
Next steps
Session information
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
##
## locale:
## [1] LC_CTYPE=pt_BR.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=es_CR.UTF-8 LC_COLLATE=pt_BR.UTF-8
## [5] LC_MONETARY=es_CR.UTF-8 LC_MESSAGES=pt_BR.UTF-8
## [7] LC_PAPER=es_CR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ohun_0.1.0 warbleR_1.1.28 NatureSounds_1.0.4 seewave_2.2.0
## [5] tuneR_1.4.1 brmsish_1.0.0 lunar_0.2-1 ggplot2_3.4.0
## [9] knitr_1.42 kableExtra_1.3.4 HDInterval_0.2.2 readxl_1.3.1
## [13] posterior_1.3.1 cowplot_1.1.1 ggdist_3.2.0 brms_2.18.0
## [17] Rcpp_1.0.9 viridis_0.6.2 viridisLite_0.4.1 remotes_2.4.2
##
## loaded via a namespace (and not attached):
## [1] backports_1.4.1 systemfonts_1.0.4 plyr_1.8.7
## [4] igraph_1.3.5 splines_4.1.0 crosstalk_1.2.0
## [7] TH.data_1.1-0 rstantools_2.2.0 inline_0.3.19
## [10] digest_0.6.31 htmltools_0.5.4 fansi_1.0.3
## [13] magrittr_2.0.3 checkmate_2.1.0 RcppParallel_5.1.5
## [16] matrixStats_0.62.0 xts_0.12.2 sandwich_3.0-1
## [19] svglite_2.1.0 prettyunits_1.1.1 colorspace_2.0-3
## [22] signal_0.7-7 rvest_1.0.3 textshaping_0.3.5
## [25] xfun_0.36 dplyr_1.0.10 RCurl_1.98-1.9
## [28] callr_3.7.3 crayon_1.5.2 jsonlite_1.8.4
## [31] lme4_1.1-27.1 survival_3.2-11 zoo_1.8-11
## [34] ape_5.6-2 glue_1.6.2 gtable_0.3.1
## [37] emmeans_1.8.1-1 webshot_0.5.4 distributional_0.3.1
## [40] pkgbuild_1.4.0 rstan_2.21.7 abind_1.4-5
## [43] scales_1.2.1 mvtnorm_1.1-3 DBI_1.1.1
## [46] miniUI_0.1.1.1 dtw_1.23-1 xtable_1.8-4
## [49] diffobj_0.3.4 proxy_0.4-27 stats4_4.1.0
## [52] StanHeaders_2.21.0-7 DT_0.26 htmlwidgets_1.5.4
## [55] httr_1.4.4 threejs_0.3.3 ellipsis_0.3.2
## [58] pkgconfig_2.0.3 loo_2.4.1.9000 farver_2.1.1
## [61] sass_0.4.5 utf8_1.2.2 labeling_0.4.2
## [64] tidyselect_1.2.0 rlang_1.0.6 reshape2_1.4.4
## [67] later_1.3.0 munsell_0.5.0 cellranger_1.1.0
## [70] tools_4.1.0 cachem_1.0.6 cli_3.6.0
## [73] generics_0.1.3 ggridges_0.5.4 evaluate_0.20
## [76] stringr_1.5.0 fastmap_1.1.0 ragg_1.1.3
## [79] oce_1.7-8 yaml_2.3.7 processx_3.8.0
## [82] pbapply_1.6-0 nlme_3.1-152 mime_0.12
## [85] projpred_2.0.2 formatR_1.11 rstanarm_2.21.3
## [88] xml2_1.3.3 compiler_4.1.0 bayesplot_1.9.0
## [91] shinythemes_1.2.0 rstudioapi_0.14 gamm4_0.2-6
## [94] tibble_3.1.8 bslib_0.4.2 stringi_1.7.12
## [97] highr_0.10 ps_1.7.2 Brobdingnag_1.2-9
## [100] lattice_0.20-44 Matrix_1.5-1 nloptr_1.2.2.2
## [103] markdown_1.3 fftw_1.0-7 shinyjs_2.1.0
## [106] tensorA_0.36.2 vctrs_0.5.2 pillar_1.8.1
## [109] lifecycle_1.0.3 gsw_1.0-6 jquerylib_0.1.4
## [112] bridgesampling_1.1-2 estimability_1.4.1 bitops_1.0-7
## [115] Sim.DiffProc_4.8 httpuv_1.6.6 R6_2.5.1
## [118] promises_1.2.0.1 gridExtra_2.3 codetools_0.2-18
## [121] boot_1.3-28 colourpicker_1.2.0 MASS_7.3-54
## [124] gtools_3.9.3 assertthat_0.2.1 rjson_0.2.21
## [127] withr_2.5.0 Deriv_4.1.3 shinystan_2.6.0
## [130] multcomp_1.4-17 mgcv_1.8-36 parallel_4.1.0
## [133] grid_4.1.0 coda_0.19-4 minqa_1.2.4
## [136] rmarkdown_2.20 shiny_1.7.3 base64enc_0.1-3
## [139] dygraphs_1.1.1.6